clc; clearvars -except NumberConsecutiveshocks LogLevel shockchoice

%  LogLevel  =  1;
%  frequency = 12; % Monthly
% %frequency =  4; % Quarterly
%-------------------------------------------------------------------------%
    load RVSP500.mat
   %sample=[ 1983 2; 2014  12];%as in  Berger, Dew-Becker and Giglio (2017), "Uncertainty shocks as second moment
   %
    sample=[ 1983 2; 2019  12];
   %sample=[ 1990 1; 2019   6]; %as in Sahalia et al. (2020) ``Volatility and Uncertainty Disconnected''
   %sample=[ 2007 7; 2009   6]; %as in Sahalia et al. (2020) ``Volatility and Uncertainty Disconnected''

% take subset
startsample=find(yyyyRV==sample(1,1)&mmRV==sample(1,2));
endsample  =find(yyyyRV==sample(2,1)&mmRV==sample(2,2));
% plot([rv*100*12 VIX])

dates_RV = yyyyRV(startsample:endsample)*100+mmRV(startsample:endsample);
RV = rv(startsample:endsample);
% Replicating Sahalia et al. (2020) ``Volatility and Uncertainty Disconnected''
fprintf('Realized volatility Market Returns %3.3f\n',mean(100*sqrt(RV*12)))

%% CHECK -- From BDBG ReStud
% load('alldata.mat','s_cme_rvar','s_cme_vixF','cme_RV','masterdates','masterdates_matlab')
% f=find(floor(masterdates/100)>=1983);
%
% load('RVSP500.mat')
% rvSP=rv;
% clear rv
% load('RVMKT.mat')
% cstr=cal(1926,7,12,size(rv));
% ss = ical(1983,1,cstr); se=ical(2014,12,cstr);
% plot([         cme_RV(f)       rvSP(ss:se)])
% plot([    sqrt(cme_RV(f)) sqrt(rvSP(ss:se))]*sqrt(12)*100)
% plot([                    sqrt(rvSP(ss:se))*sqrt(12)*100 s_cme_rvar(f)])
% plot([    sqrt(cme_RV(f))*sqrt(12)*100                   s_cme_rvar(f)])
%
% plot([log(sqrt(cme_RV)*sqrt(12)*100) ls_cme_rvar])

if LogLevel
    ls_RV = log(sqrt(RV)*sqrt(12)*100);
else
     s_RV =     sqrt(RV)*sqrt(12)*100 ;
end

clearvars -except sample NumberConsecutiveshocks RV ls_RV SP500 LogLevel dates_RV sample LogLevel shockchoice

dates_RV = dates_RV(2:end);
%% Load VIX
 VIX = xlsread('VIXCLS.xls','VIXmerged');
 VIXdates = VIX(1:end,1);
 VIX      = VIX(1:end,2); % plot([exp(ls_RV) VIX])
startsample=find(VIXdates==(sample(1,1)*100+sample(1,2)));
endsample  =find(VIXdates==(sample(2,1)*100+sample(2,2)));
 if LogLevel
     logVIX = log(VIX(startsample:endsample));
 else
     logVIX =     VIX(startsample:endsample);
 end
 clear startsample endsample VIXdates
% --------------------------------------------------------------
%% Predicting RV (Replication of )
% --------------------------------------------------------------
% % Predictive power for future 6-month RV
% 
% % Construct 6-month future var
% future_rvar=RV(2:end-5);
% for j=1:5
%     future_rvar=future_rvar+RV(2+j:end-5+j);
% end
% clear j
% future_rvar=[future_rvar;nan(6,1)];
% 
% ls_future_rvar=log(sqrt(future_rvar*2)*100);
% 
% RVBDBG_orig = ols(ls_future_rvar(1:end-6),[ones(size(ls_future_rvar(1:end-6))) ls_RV(1:end-6) logVIX(1:end-6)]);
% [[RVBDBG_orig.beta(2:end)' NaN]; [NaN NaN RVBDBG_orig.rbar]]
    
    uRVBDBG = ols(ls_RV(2:end),[ones(size(ls_RV(2:end))) ls_RV(1:end-1) logVIX(1:end-1)]);
    uRVBDBG= uRVBDBG.resid;
    uRV    = uRVBDBG/std(uRVBDBG); clear  uRVBDBG 
%corr(uRV,uRVAR12/std(uRVAR12)) = 0.9034
%corr(uRV,uRVAR1 /std(uRVAR1 )) = 0.8920
clear VIX logVIX
%% FIT AR(1)  --> Shocks 1983/M2 -- 2019/M12    
%    NumberConsecutiveshocks = 2;
% %  NumberConsecutiveshocks = 4;
% /////////////////////////////////////////////////////////////////////////
%----------- Construct indicators for seq. of shocks with 2 ident schemes-%
IndPreviousNShocksRV= zeros(size(uRV));
% /////////////////////////////////////////////////////////////////////////
for tt=NumberConsecutiveshocks:length(uRV)-1 
    if tt==NumberConsecutiveshocks
%       if (                                      sum(uEPU(tt-NumberConsecutiveshocks+1:tt)>0)==NumberConsecutiveshocks && uEPU(tt+1)<0);  IndPreviousNShocksEPU(tt) = 1; end
        if (                                      sum(uRV(tt-NumberConsecutiveshocks+1:tt)>0)==NumberConsecutiveshocks                );  IndPreviousNShocksRV(tt) = 1; end
    else
%       if (uEPU(tt-NumberConsecutiveshocks)<0 && sum(uEPU(tt-NumberConsecutiveshocks+1:tt)>0)==NumberConsecutiveshocks && uEPU(tt+1)<0);  IndPreviousNShocks2015(tt) = 1; end
        if (uRV(tt-NumberConsecutiveshocks)<0 && sum(uRV(tt-NumberConsecutiveshocks+1:tt)>0)==NumberConsecutiveshocks                );  IndPreviousNShocksRV(tt) = 1; end
    end
end
clear tt

ls_RV    =    ls_RV(2:end);

% %% To quarterly
% horizonAggreg = 3;
%   logEPUQ = logRV(13:horizonAggreg:end);
%   
%   uEPUQ = chgfreq(uRV,horizonAggreg,horizonAggreg,0); uEPUQ = uEPUQ/std(uEPUQ);
% 
% % /////////////////////////////////////////////////////////////////////////
% %----------- Construct indicators for seq. of shocks with 2 ident schemes-%
% IndPreviousNShocksRVQ= zeros(size(uEPUQ));
% % /////////////////////////////////////////////////////////////////////////
% for tt=NumberConsecutiveshocks:length(uEPUQ)-1 
%     if tt==NumberConsecutiveshocks
% %       if (                                       sum(uEPUQ(tt-NumberConsecutiveshocks+1:tt)>0)==NumberConsecutiveshocks && uEPUQ(tt+1)<0);  IndPreviousNShocksRV(tt) = 1; end
%         if (                                       sum(uEPUQ(tt-NumberConsecutiveshocks+1:tt)>0)==NumberConsecutiveshocks                 );  IndPreviousNShocksRVQ(tt) = 1; end
%     else
% %       if (uEPUQ(tt-NumberConsecutiveshocks)<0 && sum(uEPUQ(tt-NumberConsecutiveshocks+1:tt)>0)==NumberConsecutiveshocks && uEPUQ(tt+1)<0);  IndPreviousNShocksRV(tt) = 1; end
%         if (uEPUQ(tt-NumberConsecutiveshocks)<0 && sum(uEPUQ(tt-NumberConsecutiveshocks+1:tt)>0)==NumberConsecutiveshocks                 );  IndPreviousNShocksRVQ(tt) = 1; end
%     end
% end
% clear tt
% 

%#STOP
%% Pictures for paper
% [Years,M] = meshgrid(1983:2019, 1:1:12);
% calM = datenum([Years(:), M(:), 28*ones(numel(Years),1)]); % The Serial Dates
% 
% ticks = datenum([(1983:6:2019)', ones(7,1), ones(7,1) ]);
% lim_x =[datenum(  1983,01,1), datenum(2019,12,31)];
% 
% COLOR = linspecer(5); 
% LineWidth =  4.0;   
% %-------------------------------------------------------------------------%
% figure(20)
%       line(calM(3:end),   ls_RV, 'color',COLOR(1,:),'linewidth',LineWidth);      title(['log MKT volatility'])
% 
% figure(21)
% h_1 = line(calM(3:end),     uRV, 'color',COLOR(1,:),'linewidth',LineWidth);      title(['Time Series of MKT volatility shocks'])
% hold on
%       line(calM(1:end-2), 2*ones(size(calM(1:end-2),1),1),'color','k','linewidth',LineWidth);
%       line(calM(1:end-2),-2*ones(size(calM(1:end-2),1),1),'color','k','linewidth',LineWidth);     
% h_0  =line(calM(1:end-2), 0*ones(size(calM(1:end-2),1),1),'color','r','linewidth',LineWidth/2);
% uistack(h_0, 'top')
% 
% if     NumberConsecutiveshocks ==2
% figure(22)
% pos = find(IndPreviousNShocksRV==1);
% IndPreviousNShocksRV_LEAD1 = zeros(size(IndPreviousNShocksRV));  IndPreviousNShocksRV_LEAD1(pos-1) = 1;
% h_1 = stem(calM(3:end),  [uRV.*IndPreviousNShocksRV_LEAD1 uRV.*IndPreviousNShocksRV],'color', COLOR(1,:), 'LineStyle','-', 'linewidth',LineWidth/2);    
% title(['Time Series of two consecutive MKT volatility shocks'])
% %                        [uRV.*IndPreviousNShocksRV_LEAD1+uRV.*IndPreviousNShocksRV uRV]
% %     line(calM(3:end),  [uRV.*IndPreviousNShocksRV_LEAD1+uRV.*IndPreviousNShocksRV],'color', COLOR(3,:), 'LineStyle','-', 'linewidth',LineWidth)
% elseif NumberConsecutiveshocks ==3
% figure(23)
% pos = find(IndPreviousNShocksRV==1);
% IndPreviousNShocksRV_LEAD1 = zeros(size(IndPreviousNShocksRV));  IndPreviousNShocksRV_LEAD1(pos-1) = 1;
% IndPreviousNShocksRV_LEAD2 = zeros(size(IndPreviousNShocksRV));  IndPreviousNShocksRV_LEAD2(pos-2) = 1;
% [~,pos] = max(uRV); uRV(pos) = uRV(pos)/2;
% h_1 = stem(calM(3:end),[uRV.*IndPreviousNShocksRV_LEAD2 uRV.*IndPreviousNShocksRV_LEAD1 uRV.*IndPreviousNShocksRV],'color', COLOR(1,:), 'LineStyle','-', 'linewidth',LineWidth/2);    
% title(['Time Series of three consecutive MKT volatility shocks'])   
% elseif NumberConsecutiveshocks ==4
% figure(24)
% pos = find(IndPreviousNShocksRV==1);
% IndPreviousNShocksRV_LEAD1 = zeros(size(IndPreviousNShocksRV));  IndPreviousNShocksRV_LEAD1(pos-1) = 1;
% IndPreviousNShocksRV_LEAD2 = zeros(size(IndPreviousNShocksRV));  IndPreviousNShocksRV_LEAD2(pos-2) = 1;
% IndPreviousNShocksRV_LEAD3 = zeros(size(IndPreviousNShocksRV));  IndPreviousNShocksRV_LEAD3(pos-3) = 1;
% h_1 = stem(calM(3:end),[uRV.*IndPreviousNShocksRV_LEAD3 uRV.*IndPreviousNShocksRV_LEAD2 uRV.*IndPreviousNShocksRV_LEAD1 uRV.*IndPreviousNShocksRV],'color', COLOR(1,:), 'LineStyle','-', 'linewidth',LineWidth/2);    
% title(['Time Series of four consecutive MKT volatility shocks'])   
% end
% hold on
%       line(calM(1:end), 2*ones(size(calM(1:end),1),1),'color','k','linewidth',LineWidth);
% %     line(calQ(1:end),-2*ones(size(calQ(1:end),1),1),'color','k','linewidth',LineWidth);     
% hline(0)
% uistack(h_1, 'top')     
% %-------------------------------------------------------------------------%      
% 
%   set(gca, 'XLim' , lim_x)
% % set(gca, 'YLim' , [-.1 5]); %ONLY FOR FIGURES 23 & 24
% set(gca, 'XTick', ticks)
% datetick('keepticks','keeplimits')
% box on
% recessionplot
% 
% set(gca,'FontSize',32)
% %ylim([-3 6.5]) % for Figure 22
% 
% 
% %#To check consistency between AR forecasts Vs VIX+RV forecasts a-la BDBG
% IndRV_BDBG = IndPreviousNShocksRV   +IndPreviousNShocksRV_LEAD1;
% IndRV_AR12 = IndPreviousNShocksRV_AR+IndPreviousNShocksRV_LEAD1_AR;
% % They agree on sum(IndRV_AR12.*IndRV_BDBG)/2 = 26.5000  (60% of the time)
% % out of        sum(IndRV_AR12            )/2 = 49 events
% %               sum(            IndRV_BDBG)/2 = 43 events